Lmap <- function(r,X,W,cells){
  # This function calculates the vector-value of the Lorenz map of vector rank "r".
  #
  # Output of this function is used in ILF to generate a pseudo sample to calculate 
  # the ecdf, is used to calculate vector cumulative shares and to be used for the 
  # plug-in estimator for the bivariate Gini.
  #
  # Output: Lorenz map evaluated at "r"
  #
  # Inputs: r     = vector rank of interest from [0,1]x[0,1]
  #         X     = bivariate allocation (data)
  #         W     = sampling weights of X
  #         cells = the cells of power diagram pwd (pwd$cells) that signify the
  #                 vector quantile partitioning of [0,1]x[0,1]
  #
  # Forming a polygon object makes intersections easy to calculate via sf package
  N <- NROW(W)
  tryCatch({ # trycatch necessary incase cells are empty
  # rect is the bottom rectangular region [0,r_1]x[0,r_2]
  rect <- st_polygon(list(rbind(c(0,0),c(0,r[2]), c(r[1],r[2]), c(r[1],0), c(0,0)))) 
  # Check returns how the cells are intersected with rect:
  #   Case 1: Completely contained in rect
  #   Case 2: Completely outside of rect
  #   Case 3: Intersecting and not a subset
  check <- sapply(1:N, function(k){
    p <- rowProds(sweep(cells[[k]],2,r) < 0)
    return(c(all(p==1),any(p==1)&(!all(p==1))))})
  # Based on classification from "check", efficiently calculate areas of intersections
  # between cells and rect.
  return(colSums(sweep(X[check[1,],],1,W[check[1,]],"*")) + rowSums(sapply(which(check[2,]==TRUE),function(i){X[i,]*st_area(st_intersection(rect,st_polygon(list(rbind(cells[[i]],cells[[i]][1,])))))})))
  },error = function(e){})
  
  return(c(NA,NA))
}